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ABSTRACT 

We have performed cosmological N-body simulations which include the effect of the masses of the 
individual neutrino species. The simulations were aimed at studying the effect of different neutrino 
hierarchies on the matter power spectrum. Compared to the linear theory predictions, we find that 
non-linearities enhance the effect of hierarchy on the matter power spectrum at mildly non-linear 
scales. The difference between the different hierarchies is about 0.5% for a sum of neutrino masses of 
O.IeV. Albeit this is a small effect, it is potentially measurable from upcoming surveys. In combina- 
tion with neutrinoless double-/3 decay experiments, this opens up the possibility of using the sky to 
determine if neutrinos are Majorana or Dirac fermions. 
Subject headings: large-scale structure of universe — neutrinos — methods: numerical 
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1. INTRODUCTION 

Neutrino oscillation experiments have demonstrated 
that neutrinos have non-zero mass. Measurements of the 
flavor changing oscillations have provided a difference in 
the squares of the masses of the lightest and heaviest 
mass eigenstates /S.m? ~ (O.OSeV)^, yielding therefore 
a lower limit on the total neutrino mass (E = '^rn^^). 
On-going and forthcoming ground based neutrino exper- 
iments are sensitive to neutrino flavor and to the nature 
of neutrino mass (Dirac or Majorana) but are only sen- 
sitive to the absolute mass scale for large masses. On 
the other hand, cosmological probes are blind to flavor 
but sensitive to the absolute neutrino mass scale and 
there has been recently significant progress in constrain- 
ing the sum of neutrino masses from cosmological ob- 
servations. Massive neutrinos affect the observed matter 
power spectrum: their free-streaming damps the small- 
scale power and modifies the shape of the matter power 
spectrum below the free-streaming length ( Doroshkevich 
' —^^^ — — — ~ - - - ^^2Q06j ^Kiako- 



et al.||I980t |Hu et al.||1998l [TakadiTetal 



tou et al.|2008p . Up-to-date only upper limits have been 
obtained. However, the constraints are getting tighter 
-the c urrent limits o n the total mass are < 0.3eV (e.g., 
[Thomas ct al. (2009); Reid ct al. (2010a); Komatsu et al. 



(|20IfJ;.Saitoetal..(.20fI 
deTIitteret al.J ( |2012p )- 



j; ,Riemer S0rcnscn ct ar "(2011 i 



and closer to the experimental 
limits derived from accelerator, reactor, solar, an d at- 
mospheric neutrino oscilla ti ons (see the reviews by Les- 
purgues fc Pastor] (2006); Gonzalez-Garcia fc Maltom 

Detecting the effect of 
neutrino masses on cosmological structure and measur- 
ing the neutrino mass scale is well within the re ach of up- 
coming cosmologica l surve ys 



•:S^ 



Takadaet al. (2006); 



iHannestad & Wong' (2007'); 'Kitching et al.'f^OOS); LSST 
1|009); Hanncs tad (2010); Lahav ct al. (2010); Rcid et aT 
^OlOa); Juncn ez et al.| l |20lO| ); |Carbone et aT] i |2011[ ) anH^ 
references therein) . 

The neutrino mass splitting required to explain the 
oscillation results implies that for three neutrino species 
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there are two possible hierarchies in the neutrino mass 
spectrum: normal hierarchy (NH) with two light states 
and one heavy one and a total mass E > 0.05eV; inverted 
hierarchy (IH) with two heavy states and one light one 
with E > O.IeV. If the absolute mass scale is much higher 
than the mass splitting then the mass hierarchy does not 
matter and this case is referred to as degenerate mass 
spectrum. The degenerate hierarchy however is under 
pressure from observations (Thomas et al. 2009[ Rejd 
et al. 2010a Gonzalez-Garcia et al. '2010); see Fig. [T] 



where we have introduced the neutrino mass splitting 
parameter A relating the lightest (m) and heav iest [Ad] 



neutrino masses following Jimenez et al. ( 20f0| ) [^ 

NH : E = 2m -h M A = (M - m)/E 
IH: E = m + 2il/ A = (m-M)/E. 



(1) 



To determine the possible hierarchy is very relevant 
as it can complement results from neutrinoless double-/? 
decay to help determ ine the nature of the neutrino itself 
(Jimenez et al. 2010): is it its own anti-particl e? that 
is, IS i t "a~M aiorana fermion? As discussed in |Bahcall| 
et al. (2004) (see their Fig. 3 recalling that a light neu- 
trmo mass of 0.07eV corresponds to E ~ 0.2 — 0.25eV), 
if the next generation of neutrinoless double-/3 decay ex- 
periments find a signal, then neutrinos are Majorana. If 
these experiments do not see a signal, it is important to 
discriminate if that is because the signal is below the de- 
tection threshold, or neutrinos are truly Dirac particles. 



Here is where cosmology enters (Jimenez et al. 2010): if 
E > 0.25eV, then the hierarchy is eftectively degenerate 
and neutrinos are Dirac. The interesting region to de- 
termine the hierarchy is for O.leV < E < 0.25eV, where 
the absence of neutrinoless double-/? decay indicates that 
neutrinos are Dirac o nly if the hierarchy is inverted. In 
Jimenez et al. (2010) we addressed the above question 
using linear tlieory to predict the effect of neutrinos on 
the matter power spectrum. Our conclusion was that an 



ambitious survey a la stage IV dark energy task force ( Al 
brecht et al.|2006 1 , could in principle distinguish between 
the invertecl and normal hierarchy if E < 0.1 5eV. How- 

■^ Since one mass splitting is much larger than the other, for cos- 
mological applications we can safely ignore the small mass splitting. 



ever, because the effect on the matter power spectrum 
of neutrinos extends into the non-hnear regime, it is cru- 
cial to perform N-body simulations that include massive 
neutrinos to properly quantify the effect. 

The main physical effect that distinguishes differ- 
ent hierarchies is the fact that neutrinos of different 
masses have different transition from relativistic to non- 
relativistic thus influencing the shape of the matte r 



power spectrum ( Slosar|[2006 de Bernardis et al. 2009). 
Note that neutrino oscillatioii experiments rule out large 
regions in the E — A parameter space (see Fig.jll) and 
therefore it is worth investigating only the allowecfregion 
(gray swath in Fig. [I]) . 
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Figure 1. Constraints on the mass splitting from neutrino oscil- 
lations (shaded regions) and total neutrino mass from cosmology 
(vertical dashed line) in the parameter space defined by the sum 
of neutrino masses S and the mass splitting parameter A charac- 
terizing the hierarchy. The key region where it is interesting to 
determine A is 0.1 eV < S < 0.3 cV. The triangles indicate the 
cases studied in this paper for normal and inverted hierarchies (at 
S = O.leV) and for th e degenerate hierarchy (at S = 0.3eV). Plot 
adapted from Fig. 1 of|Carbone et al.]l|2011[|. 



2. NUMERICAL SIMULATION METHOD 

In this Letter we study the effect of the neutrino mass 
splitting on the non-linear matter power spectrum. In 
particular, we are interested in the question of how and 
if the neutrino hierarchy leaves an imprint on the matter 
power spectrum in the non-linear regime. To address this 
question we run cosmological N-body simulations for sev- 
eral different neutrino mass configurations: A) 3 massless 
neutrinos; B) 3 neutrinos of equal mass (degenerate case, 
E = 0.3eV); C) 1 massless neutrino and 2 massive neu- 
trinos (inverted hierarchy, E = O.leV, A = —0.50); D) 2 
light neutrinos and 1 heavy neutrino (normal hierarchy, 
E = O.leV, A = 0.32). 

So far there have been mainly two approaches to incor- 
porating neutrinos into cosmological N-body simulations: 
sampling the neutrino densit y with particles just like in 
the case of cold dark matter (White et al.' 1983; Klypin 
et al.||1993| iBrandbyge et al. 2008; Viel et al. 20l^ 



JC 



or 



evolving the neutrin o density on a grid (w ith a hxcd~ size) 
using linea r theory ( Brandbyge fc IIannestad||2009 Viel 
et al.|2010 ). Thes e two approaches were cornbmed into a 
hybrid method by Brandbyge & Hannestad (2010). The 
particle-based approach has the advantage or beiiig able 



to capture the non-linear evolution of the neutrinos and 
their effect on the cold matter components, while the 
linear evolution of the neutrino density on a grid is by 
construction only able to model the linear gravitational 
effect of the neutrinos on the non-linear matter distribu- 
tion. In the particle-based approach, however, the neutri- 
nos are always treated as non-relativistic particles, since 
standard cosmological N-body codes do not include rela- 
tivistic effects. Another problem of this approach is that 
the finite number of particles used to sample the neutrino 
density generate shot noise. As the neutrinos have high 
thermal velocities, they move quickly away from their ini- 
tial positions and give rise to a Poisson-like shot noise. 
Both these problems can be somewhat circumvented by 
starting the simulation at a sufficiently low redshift, when 
the neutrinos are practically non-relativistic and their in- 
put power spectrum is large compared to the shot noise. 
On the other hand, one does not want to start too late, 
when non-linearities have already become significant on 
the scales of interest. Hence, in order to keep the shot 
noise sub-dominant already at the initial redshift of the 
simulation (ideally this would be the redshift at which 
the neutrinos become effectively non-relativistic), a large 
number of n eutrino tracers is required 



S tudies by Brandbyge fc Hannestad] ( 2010 ) and |Bird"et 
al. (2012) have shown that in spite oif the shortcomings 



mentioned above the particle-based approach is more ac- 
curate in computing the effect of massive neutrinos on the 
non-linear matter power spectrum than the grid-based 
approach. This is especially the case, when looking at 
relative quantities like the ratio of the power spectrum 
with and without massive neutrinos, since in this case 
systematic effects due to the rather late start of the sim- 
ulation cancel out to a large extent. Hence, in this paper, 
we present results of particle-based simulations only. 
We focus on scales in the range of 0.01 /iMpc~ < k < 

1 hMpc~^ . These scales will be probed to high accuracy 
by future surveys and are much less affected by baryonic 
physics (which we neglect in this paper) than smaller 

scales. Hence, we simulate a volume of (600/i~^Mpc) 
with a particle load of 1 billion cold matter tracers and 

2 billion neutrino tracers. 

We adopt a flat ACDM cosmology with cosmological 
parameters compatible with current observational con- 
straints. The primordial curvature power spectrum is 
specified by the scalar amplitude A^ — 2.45 x 10^^ at 
the pivot scale kp = 0.002 Mpc~^ with a spectral index 
Us = 0.97. We keep the present-day total matter frac- 
tion the same for all neutrino models: flm = ^cdm + 
fl^ + rit = 0.27 with the baryon fraction and the neutrino 
fraction given by ^b = 0.046 and n^ = E/(93.8eV/i2), 
where E is the sum of neutrino masses in units of eV. 
We choose the Hubble parameter to he h — 0.7. This 
choice of cosmological parameters yields a present-day 
linear mass variance in spheres of 8 h^^ Mpc of trg ~ 0.8. 

For setting up the initial conditions and for assessing 
the N-body simulation results, we need accurate linear 
predictions for the transfer functions for each compo- 
nent of the universe. To this end, we use the linear 
Einstein-Boltzman n solver CAMB (Version Jan 2012) 
( Lewis et al.|2000" |, which computes the individual linear 
transfer functions for cold dark matter, baryons, neu- 
trinos, and photons. As we simulate the different neu- 
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Figure 2. Fractional difference in the total matter power spectrum (CDM+baryons+massive neutrinos) of the massive neutrino runs 
and the massless neutrino run. Data points show the simulation results. The solid lines depict the linear theory predict ions, while the 
black dotted lines show the estimated effect including non-linearities using the extended HALOFIT formula (Bird et al. 2012|. The left 
panel shows the degenerate case with S = 0.3eV. The middle and right panel show the normal and inverted hierarchy, respectively, with 
S = O.leV. Note the typical shape of the effect in the non-linear case: non-linearities enhance the effect and lead to a maximum in the 
suppression located at mildly non-linear scales. 



trino species separately with different neutrino tracers, 
we modified CAMB to store for each neutrino species a 
separate transfer function. 

To reduce transients d ue to the late start of the simula- 
tion ([ Scoccimarro|1998 ), we implement second-order La- 
grangian perturbation theory (2LPT) for the cold matter 
components (cold dark matter and baryons). We com- 
pute numerically the (fc-dependent) linear growth rate 
from two CAMB outputs around the initial redshift. 
We approximate the second-order growth function D2 
by D2 ~ —3/7D^ w here Di is linear growth function 
( Bouchet et al.|1995 ). At the initial redshift (z^ = 9) the 
neutrino perturbations are still very much in the linear 
regime. Hence, for the neutrino particles, it suffices to 
use the Zel'dovich approximation to displace the parti- 
cles from the initial grid points and to assign the gravi- 
tational velocities. Then thermal velocities drawn from 
the appropriate Fermi-Dirac distribution are added to 
the neutrino velocities. We neglect higher-order multi- 



poles of the neutrin o phase space distribution (see Ma & 
Bertschinger ( 19941) for a treatment of the full neutrmo 
phase space), i'ests with CAMB where these multipoles 
were set to zero at z = Zi, and were then evolved fur- 
ther to low redshift, have shown that the effect of these 
initial conditions o n the linear matt er power spectrum is 
negligible at 2 < 2 ( |de Putter||2012[ ). 
The simulations were carried out with Gadget- 2 



( [Spring cl 2005). We modified Gadget-2 shghtly to take 
into account the effect of the massive neutrinos (and the 
radiation) on the evolution of the scale factor a. 

Using smaller test runs, we performed convergence 
tests on the neutrino time stepping. We found that a 
maximum time step max(Alna) — 0.025 in the long- 
range particle-mesh force computation is sufficient for 
an accurate computation of the effect of the neutrinos on 
the matter power spectrum. Note that we disabled the 
Courant condition for the neutrino tracers. Hence, the 
time step for the long-range force is determined from the 
velocity dispersion of the cold matter tracers alone. This 
speeds up the computations significantly and leaves the 
non-linear matter power spectrum virtually unchanged. 
We set the softening length of the short-range force to 
20 kpc/h for both the cold matter and neutrino compo- 
nent. 

Varying the initial redshift and the number of neutrino 



tracers, we confirmed that the measured ratio of the non- 
linear matter power spectrum with and without massive 
neutrinos is robust against the neutrino sho t noise and 
the residua l transients due to the late start (Brandbyge 
eraL][2008| ) 



We compute the total matter density contrast S by as- 
signing the particles to a 2048'^ grid using the cloud-in- 
cell scheme. In this process cold matter and neutrino par- 
ticles are weighted by the fraction they contribute to 51™. 
Using Fast Fourier transforms and averaging \5k\'^ over 
spherical shells with a bin width of Afc — 0.01 h/Mpc, we 
obtain the power spectrum P{k) = (|(5fcp). Note that we 
only consider the non-relativistic neutrino species sam- 
pled by particles, i.e. the fluctuations in the radiation 
(relativistic neutrinos and photons) are not taken into 
account. At the redshifts of interest (z < 2), however, 
radiation contributes a negligible amount to the total en- 
ergy budget. One possible uncertainty in our simulation 
results is the shot noise coming from the cold matter par- 
ticles. As the particles start off from a grid, the shot noise 
is sub-Poissonian and scale-dependent at high redshifts. 
However, for z < 1 and k < lft,Mpc~^, even a Poisson 
shot noise is negligible, 
correct for shot noise. 



Hence, we do not attempt to 



3. RESULTS 

First, we compare the results of the degenerate case 
with S — 0.3eV with previous works. We consider the 
fractional difference in the total matter power spectrum 
of models with and without massive neutrinos. The ad- 
vantage of this relative quantity is that the sample vari- 
ance present in the N-body simulations cancels out al- 
most completely. The suppression of the total matter 
power spectrum due to the three degenerate massive neu- 
trinos is shown in the left panel of Fig. J2j The simula- 
tion results show that non-linearities enhance the effect 
on mildly non-linear scal es. This behavior is anticipated 



by perturbation theory (Saito et al. 2008 Wong 2008 



[Lesgourgu es et al.|20 09) . In contrast to the linear theory 
predictions (solid lines), in the non-linear case there is 
a maximum suppression, whose depth and position in k 
depend on redshift. These numerical results are in ex- 



cellent agree ment with Brandbyge et al. ( 2008 ) and Viel 
et al. (12010 ) and can be reasonably well modeled with 
HALOFrri lSmith et aL|,2003) , which was extended by 




0.1 
k [h/Mpc] 



Figure 3. Fractional difference in the total matter power spec- 
trum (CDM+baryons+massive neutrinos) of the inverted hierar- 
chy and normal hierarchy run (the sum of neutrino masses is kept 
fixed and only the mass splitting is varied). Note that also in this 
case, non-linearities enhance the effect on mildly non-linear scales 
compared to linear theory predictions (solid lines). 



Bird et al. ( 2012 ) to model the effect of massive neutrinos 
(black dotted lines) . 

In the middle and right panel of Fig. [2] we show the 
power spectrum suppression for the normal and inverted 
hierarchy with S — O.leV. We observe the same qual- 
itative behavior as before. In order to make the small 
difference between the two models visible, the fractional 
difference in the total matter power spectrum of the two 
cases is shown in Fig. [3] Similar to the case of mas- 
sive vs. massless neutrinos, the non-linear evolution en- 
hances the difference between the two models on mildly 
non-linear scales. On even smaller scales which are in 
the stable clustering regime, the strong non-linearities 
eventually overcome the initial differences between the 
models and the remaining effect decreases with fc and 
may eventually drop below the linear theory prediction 
at large fc. 

In order to test if the simulation results are affected by 
the realization of the initial Gaussian random field, we 
compare them with the results of another set of simula- 
tions with a different random realization. We find that 
both sets of simulation give virtually the same predic- 
tions. 

For completeness, we show the neutrino power spectra 
for the normal and inverted hierarchy in Fig. |4] On large 
scales (A: < 0.1 /i/Mpc), the neutrino power spectra from 
the simulations agree well with the linear predictions. As 
expected, even at high redshift we observe a Poisson-like 
shot noise due to the large thermal velocities of the neu- 
trinos. Although this shot noise is large in the neutrino 
power spectrum, it is much smaller when one considers 
the total matter density, since neutrinos contribute less 
than 1% to the total matter fraction. Additionally, on 
scales in the non-linear regime, the matter power spec- 
trum is dominated by the non-linear power sourced from 
large-scale modes, for which shot noise is negligible. 

4. DISCUSSION AND CONCLUSIONS 

The main result of our numerical experiments is that 
non-linearities enhance the dependence of the power 



spectrum on the different neutrino hierarchies, thus mak- 
ing the observational signature more pronounced. We 
estimate that, if all other cosmological parameters are 
known (including the sum of neutrino masses E), the 
two hierarchies can be distinguished with confidence, as 
illustrated in Fig. [5] as function of the maximum k con- 
sidered, making the effect potentially measurable. We 
have assumed an effective volume of 1 (Gpc//i)'^ at z = 
(red lines) and 10 (Gpc//i)^ at z = 1 (blue Hnes)Q 

Whether degeneracies with other cosmological parame- 
ters and systematic effects (galaxy bias, baryonic physics, 
observational limitations etc.) will cancel the detectabil- 
ity of the effect re mains to be explored and will be con- 
sidered elsewhere ( Wagner et al.|[20T2 ). 

Our findings indicate that cosmology has the potential 
of determining the neutrino hierarchy in the interesting 
window S > O.leV. Signal-to-noise estimates done using 
linear theory predictions indicated that if S happens to 
be in the window O.lSeV < S < 0.25eV, cosmology could 
not help determi ne the hiera rchy and thus the nature of 
neutrino masses (Jimenez et al. 2010), leaving an impor- 
tant gap in our knowledge of neutrino properties. The 
fact that non-linearities enhance the effect compared to 
the linear prediction, will potentially enable cosmology 
to determine the hierarchy in a wider E range and pos- 
sibly close the gap. 



As an aside and already noted in Jimenez et al. (2010), 
cosmology is more sensitive to |A| than to its sign: a 
measurement of |A| in agreement with that predicted by 
oscillations experiments for the measured E would pro- 
vide a convincing consistency check for the total neutrino 
mass constraint from cosmology. 

However, it is very important to keep in mind that one 
needs theoretical predictions of the absolute non-linear 
power spectrum at least accurate to the 0.1% level to ac- 
tually be able to make any sensible measurement of the 
neutrino hierarchy. In this Letter we presented numerical 
predictions only for the relative non-linear power spec- 
trum suppression. This relative quantity is much more 
robust against numerical errors. Even without massive 
neutrinos, it is challenging to compute the non-linear 
power spect rum to sub-percent precision (e.g.,|IIeitmann 



et al. 20101). On small scales {k > 1/i/Mpc), baryon 
physics, which is strongly model dependent and com- 
putationally very intensive, is non-negligible and makes 
it very hard to achieve this accuracy in the foreseeable 
future. In addition, although it is in principle much 
less demanding to compute the linear power spectrum to 
high accurac y, different linear E instein-Boltz mann codes 
(e-g. CAMB dLewis et al.|[200iol ) and CLASS ( Bias et al. 
2011 )) still do not agree to 0.1% precision on the relevant 
scales. 

Despite these very challenging and open problems, 
precision measurements of the large-scale structure of 
the universe remain an interesting avenue to determine 
the neutrino hierarchy. 

We acknowledge the participation of Carlos Pcha- 
Garay at an early stage of this project and thank him 
for many interesting discussions. CW is grateful for in- 

* These volumes roughly correspond to the volume out to ^ = 0.5 
and between z = 0.5 and z = 1.5 in 1/10 of the sky respectively in 
a standard ACDM universe. 
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Figure 4. Neutrino power spectra for the normal (left panel: light neutrino, middle panel: heavy neutrino) and inverted hierarchy (right 
panel) at several redshifts. The shot noise effect described in the text can be easily seen. In our approximation and for the chosen value 
of E, the normal hierarchy has effectively two non-zero neutrino masses (M and m in the notation of Eq. 1), while the inverted hierarchy 
has effectively only one massive eigenstate M (the light neutrino is massless). 
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